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ABSTRACT 

The universality of the halo mass function is investigated in the context of dark energy cos- 
mologies. This widely used approximation assumes that the mass function can be expressed 
as a function of the matter density and the root-mean-square linear density fluctuation 
a only, with no explicit dependence on the properties of dark energy or redshift. In order to 
test this hypothesis we run a series of 15 high-resolution N-body simulations for different 
cosmological models. These consist of three ACDM cosmologies best fitting WMAP-1, 3 
and 5 years data, which are used for model comparison, and three toy-models characterized 
by a Ratra-Peebles quintessence potential with different slopes and amounts of dark energy 
density. These toy models have very different evolutionary histories at the background and 
linear level, but share the same cts value. For each of these models we measure the mass func- 
tion from catalogues of halos identified in the simulations using the Friend-of-Friend (FoF) 
algorithm. We find redshift-dependent deviations from a universal behaviour, well above nu- 
merical uncertainties and of non-stochastic origin, which are correlated with the linear growth 
factor of the investigated cosmologies. Using the spherical collapse as guidance, we show that 
such deviations are caused by the cosmology dependence of the non-linear collapse and viri- 
alization process. For practical applications, we provide a fitting formula of the mass function 
accurate to 5 percents over the all range of investigated cosmologies. We also derive an empir- 
ical relation between the FoF linking parameter and the virial overdensity which can account 
for most of the deviations from an exact universal behavior Overall these results suggest that 
measurements of the halo mass function at z = can provide additional constraints on dark 
energy since it carries a fossil record of the past cosmic evolution. 

Key words: N-body simulations, dark energy, mass functions, large-scale structures, 
quintessence, cosmology 



1 INTRODUCTION 

In this series of articles we have investigated the imprint of dark 
ener gy on the non-lin ear structure formation. In a previous pa- 
per ( lAlimi et al. I2OIOI) we focused on the non-linear matter power 
spectrum at z = and showed that dark energy leaves distinctive 
signatures through a number of effects. On the one hand the cluster- 
ing of dark energy modifies the shape and amplitude of the linear 
matter power spectrum, on the other hand the values of the cos- 
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mological parameters, such as the matter density $1™, and the lin- 
ear root-mean-square (rms) density fluctuations on the 8 h^^Mpc 
scale, (T8, differ from one model to another such as to satisfy the 
constraints from Supernova la and Cosmic Microwave Background 
(CMB) observations. Because of this, the shape and amplitude as 
well as the evolution of the linear power spectrum are affected, with 
the non-linear phase of collapse mixing and amplifying these model 
dependent features. This is a direct consequence of the fact that a 
record of the past history of forming structures is kept through- 
out the non-linear regime (see e.g. Ma 2007). Although current 
measurements of the clustering of matter at small scales are unable 
to detect such imprints (mainly because of astrophysical system- 
atic uncertainties related to galaxy bias), these effects are present 
and may become detectable with future weak lensing observations 
(see e.g. ?). Another consequence is that any estimate of the non- 



2 J. Courtin, Y. Rasera, J. -M. Alimi, P. -S. Corasaniti, V. Boucher, A. Fuzfa 



linear power spectrum based on parametrized fitting functions writ- 
ten in terms of the cosmological parameters (e.g. Qm) and instan- 
taneous linear quantities (such as the linear matter power spectrum, 
P{k), at z = and the linear growth factor) are of limited pre- 
cision on the non-l i near scales. For ins t ance t his is the case of the 
ISmith et al. 1 ( l2003h : |Peacock & Doddj ( Il996h formula. Deviations 
with respect to these fitting functions depend on the past evolution- 
ary history of a given cosmology, hence it is not surprising that 
such discrepancies have been found to be manifestly accentuated 
in the context of dar k energy cosmologies I McDonald et al. IllOOd: 



Ma'2007':'Fra ncis et al. l2007l : ICasarini et al. l2009l : I Jennings et all 
2010; Alimi et al. Il2010h . 

The halo density profile i s another observable fo r which sim- 
ilar effects occur. For example IWechsler et al. I ^2003) have shown 
that the concentration parameters depend on the halo assembly his- 
tory, and it has been shown that such dependencies are strengthened 
in the case of dark energy models (see for instance Dolag et al. 
l2004h . Paradoxicrily, as a result of state-of-the-a rt numerical simu- 
lations jjenkins et al7ll200ll : lmrren et al. II2006I) . a universal form 
of the halo mass function entirely specified by Sl^ and linear root- 
mean-square (rms) density fluctuations a is usually assumed. Fur- 
thermore it has been claimed that such a universal form h olds for 
dark energy cosmologies as well dLinder & Jenkins1l2003h . If uni- 
versality is rigorously exact, then in the light of the previous results 
on the non-linear matter power spectrum and halo profile, it would 
imply that there must exist an unknown gravitational mechanism 
capable of erasing the influence of the past evolution of forming 
structures on the halo mass function. Only in such a case a de- 
pendence on cosmology (e.g. the properties of dark energy) and 
redshift would be absent. 

The universality of the mass function at very high-redshift 
(z > 5) has been long debated in the literature. As an exam- 
ple devia tions from a universal be haviour up to 50% have been 
found by iReed et all j2003l.l2007l) . However iLukic et al. I ilOOH) 
have shown that most of these deviations might be caused by nu- 
merical artifacts (finite volume effects or initial conditions set at 
very low redshift). In fact results inferred from high redshift simu- 
lations are very sensitive to numerical errors, consequently several 
wo rks have focused on t he low-redshift mass function. As an exam- 
ple [Tinker etaTI ( l2008h have shown deviations from universality 
up to 30% in the low redshift mass function (z < 3). Such devi- 
ations are thought to be associated with the Spherical Overdensity 
(SO) halo finder which tends to underestimate the mass function 
relative to the F riend-of-F riend (FoF) algorithm when considering 
higher redshifts jLukic et al. 20 09). Since the SO detection is sim- 
ilar to th e observational proc edure of measuring the mass of galaxy 
clusters jTinker et aini2008 l). such an effect is more relevant from 
an observational stand point, but less informative for a better un- 
derstanding of the non-linear structure formation. To our knowl- 
edge, 10% deviations from universality (as a function of redshift ) 
using FoF halo find er have been detected by iLukic et al. 1 ( l2007h : 



Tinker et al. at low redshift and very recently confirmed by 

Crocce et al. (2010). Nonetheless the physical origin of these de- 
viations has yet to be understood, especially in the context of dark 
energy cosmologies. 

Is the halo mass function really universal? To what extent does 
the universality approximation hold? For which cosmologies and 
for which redshifts? If there are deviations from an universal be- 
haviour are these of stochastic nature or do they correlate with 
physical effects? What can we learn from deviations to a universal 
behaviour and how to model them? These are the questions which 
we will address in this work. 



The paper is organized as follows. In Section|2]we introduce 
the halo mass function, describe the main features of the cosmo- 
logical models for which we have run a series of N-body simula- 
tions, and discuss the spherical collapse model. In Section [3] we 
describe the characteristics of the N-body simulations, and the halo 
finder algorithm, while in Section |4l we discuss various numeri- 
cal tests which we have performed to identify potential sources 
of systematics errors. In Section [5] we present the results of the 
non-universality of the mass function, and discuss the mechanisms 
responsible for the measured deviations in Section |6] We finally 
discuss our conclusions in Section|7] 



2 DARK ENERGY AND STRUCTURE FORMATION 
2.1 The halo mass function 

Current analytical predictions of the halo mass function are based 
on the original work by Press & Schechter ( 1974). The basic idea 
is that virialized objects of mass M correspond to regions where the 
linear density fluctuation field smoothed on the scale R lies above 
a critical density contrast threshold Sc. Then the halo mass function 
is simply proportional to the fraction of volume occupied by the 
collapsed objects with mass greater than Al. Assuming a Gaussian 
distribution of density fluctuations, this is given by: 



where 



F{> M) = 



27rcr 



k^P{k)W'^{k,R)dk, 



(1) 



(2) 



is the variance of the density field smoothed on the scale R, 
with P{k) being the linear matter power spectrum today and 
W{k,R) is the window function. For a spherical top-hat filter 
in real space of radius R containing a mass M ~ Aivpo/SR^ 
where po is the present matter density, we have W{k,R) = 
3 X {sin{kR) I (kRf - cos{kR)/{kRf). The only additional in- 
gredient needed to solve the model is the overdensity threshold 5c, 
which is assumed to be given by the spherical collapse model pre- 
diction of the linearly extrapolated density fluctuation at the time 
of collapse. 

Subsequent studies of the mass function h a ve largely im- 
prove d this simple modelling l lBond et al. I Il99ll : iLacev & Cole I 
Il993h . and corrections have been included to account for the 
ellipsoidal collapse lAudit et al. 1 1 19971 : ISheth & Tormen 1 1 19991: 
Sheth et a l. 2001). Recently a modelling of the halo mass function 
in the context of the excursion set formalism as well as its extension 
to the case of non-gaussian initial conditions has been presented in 
jMaggiore & Riotto II2009I) . 

A clear understanding of what universality of the mass func- 
tion implies can be gathered by writing Eq. lO as 



F{> M) 



d 



(3) 



where S{5,S/cr) is a selection function in 5-space and is 
the probability distribution of the primordial density fluctuations 
smoothed o n the scale R (for more details about the following dis- 
cussion see iBlanchard et al7lll992l) . In such a case the halo mass 
function reads 



dn 
dM 



1 da 



M cr2 dM 



d5. 



(4) 
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This very general formulation should be a good approximation for 
slowly evolving primordial power spectra (such as the one used in 
this paper). A crucial point resulting from Eq. ((U is that all effects 
associated with the non-linear collapse are encoded in the form of 
the selection function S. Although the precise shape of S is difficult 
to compute, one may expect as a general trend that S varies from 
zero to one near the non-linear density threshold Ss- In the end it is 
this threshold that determines the precise form of the mass function. 
For instance, by assuming 5 to be a Heaviside in 5s, one recovers 
the Press-Schechter halo mass function. 



dn _ po 1 d(j f 5s\ )■ 



(5) 



The dependence on Ss accounts for the effect of the non-linear 
gravitational collapse and the virialization process. The former is 
estimated in terms of the extrapolated linear density at the time of 
collapse 5^, and the latter by the virial overdensity at the same time 
Afi Jj. Since the collapse and virialization processes are specific to 
each mass, redshift and cosmology, we can expect 5s to be cos- 
mology and redshift dependent, and consequently, for a given a 
the mass function as well. This is explicitly manifest in the Press- 
Schechter (PS) mass function formula, with 5s = <5c 



dn p 1 da ( 5c 

Im " 'JI'^dM-' W 



where the functional form of / is given by: 



fps 



5c 



1/2 



(6) 



(7) 



While, in the c ase of the Sheth - Torme n 
jSheth & Torriien] 1 19991 ; Isheth et al. I 1200 ll) . 
form of / is given by: 



formula (ST) 
the functional 



fsT 



5c 



A 



2a 



1/2 



1 + 



5c 



-2p 



(8) 

with A = 0.322, a = 0.707 and p = 0.3. Again in these 
formula 5c is the spherical collapse prediction of a given cos- 
mology, nevertheless the cosmology dependen ce encoded in 5 c 
is often negle c ted in the literature (see how ever Percival 1 ( l2005h : 
iFrancis et al. I j2009al lbl): iGrossi & Springel ( l2009h ). This is be- 
cause the spherical collapse model in the SCDM scenario, which 
has been for long time the reference cosmology to study structure 
formation, predicts 5c = 1.686 and A^^^ ~ 178 constant in red- 
shift. Therefore it has become common to set 5c to such a value. Al- 
ternatively the functional form of / in Eq. ^ has been directly fit- 
ted against the mass function measured in numerical simulations as 
a function of cM3nly (Jenkins et al. 2001; Linder & Jenkins 2003; 
IWarren et al7ll2006h . For example, I Jenkins et al. I j200lh found 



/(cr) = 0.315 exp(-| Incr 



+0.611^-*) 



(9) 



over the range —1.2 < Incr"^ < 1.05, with deviations for differ- 
ent cosmologies within 20% level. In such a case Eq. ^ depends 
on the matter density Qm and the rms linear density fluctuation a 
only, thus manifestly independent of the specificities related to the 
cosmological and redshift evolution of the non-linear collapse and 
virialization process. This is what is commonly understood as ''uni- 
versality'' of the mass function. Another way of rephrasing this idea 



^ These quantities do not refer to any specific non-linear collapse model, 
thus we distinguish them from 5c and Ayir of the spherical collapse. 
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Figure 1. Relative variation of f{5c/o-) for different values of the spherical 
collapse derived density threshold Sc as a function of ln((T^^). Top panel 
shows variation for the Press-Schechter formula, while we plot the varia- 
tion for the Sheth-Tormen formula in the bottom panel. Curves from top to 
bottom con'espond to regularly spaced values of 5c in the range 1.638 to 
1.686 (covering the interval of values inferred from the spherical collapse 
for the cosmological models considered in this paper). 



is to say that the function / in Eq. ^ 



's^H-)dS, 

o5 a 



(10) 



is universal if the selection function S is independent of cosmol- 
ogy and redshift. How exact is this statement, and to what extent it 
remains valid especially in the context of dark energy cosmologies? 

We can have an estimate of the influence of varying 5s on 
/((j) by evaluating the relative error 



f 







as- 



'•-Hi)d5 



rex. 

Jo 



I d5 



A5s, 



(11) 



we may notice that deviations from a universal behaviour are 
proportional to A5s. In Fig. [T] we illustrate this for the (ex- 
tended) Press-Schechter (PS) formula and the Sheth-Tormen (ST) 
parametrization, where we let 5s (= 5c) varying in the range 1.638 
to 1.686 (corresponding to the typical range of 5c values covered 
by the cosmological models studied in this paper). We can see that 
the effect of varying 5s is more important on the high mass tail of 
the mass function (i.e. ln{a~^) > 0). This is because the mass 
function is exponentially sensitive to the value of 5s ■ It is for this 
very reason that we will specifically focus on the high mass tail, 
nevertheless this does not imply that there are no effects on smaller 
masses. From Fig.[T]we can also notice that the variations induced 
on f{(j) are similar for the PS and ST prescriptions, thus suggest- 
ing that we can study deviations from universality independently 
of the specific form fo the mass function. To this purpose, we focus 
on models which exhibit the same rms fluctuations a, while be- 
ing characterized by different background expansion histories and 
evolutions of the linear perturbations. 
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Model 


^DE 


h 


0-8 




Hi, 


ACDM-W5 


0.74 


0.72 


0.79 


0.963 


0.044 


ACDM-W3 


0.76 


0.73 


0.74 


0.951 


0.042 


ACDM-Wl 


0.71 


0.72 


0.90 


0.99 


0.047 



Table 1. Cosmological parameter values of the WMAP calibrated cosmolo- 
gies. Our reference model is the ACDM-W5. 

2.2 Cosmological toy-models: background and linear 
evolution 

We consider as refer ence cosmology a A CDM model best fit to 
WMAP-5 years data jKomatsu et al. 112009:) which hereafter we re- 
fer to as ACDM-W5. For comparison we also consider two ad- 
ditional ACDM models cal ibrated to WMAP-1 and 3 years data 
dSpergel et afl |2003|. |2007|) . which we dub as ACDM-W3 and 
ACDM-Wl respectively. The model parameter values of these 
ACDM- WMAP cosmologies are listed in Table [T] the most no- 
ticeable difference between these models concerns the erg value, 
nevertheless their expansion and linear growth histories are almost 
identical as we will discuss later in this Section. In order to inves- 
tigate the imprint of dark energy on the halo mass function and 
test the universality hypothesis, we confront the ACDM-W5 cos- 
mology with a set of "toy-models". These are flat cosmological 
models with different background expansion and linear growth of 
the density perturbations. Following the discussion of the previous 
paragraph we additionally require the models to have the same dis- 
tribution of linear density fluctuations at z = 0, hence the same as 
value. 

We focus on a q uintessence model with Ratra-Peebles poten- 
tial ( lRatra& Peebles ..1988 ): 

^W) = ^7-7^' (12) 

where a and A are the slope and amplitude of the scalar self- 
interaction respectively, and (fi is the quintessence field evolving ac- 
cording to the Klein-Gordon equation. For a = the quintessence 
cosmology resembles a standard ACDM, provided that the initial 
field velocity vanishes. We choose the RP model since it corre- 
sponds to a dark energy component whose equation of state can 
vary from w = —1 (cosmological constant value) to an evolving 
function of the redshift (w{z) > —1) by changing the slope of the 
potential through the parameter a > 0. As we are specifically in- 
terested in cosmological evolutions which largely differs from that 
of ACDM, we consider a quintessence model with a = 10, which 
we dub as L-RPCDM (the letter L in the acronym means large a 
value). 

We also construct two models with different amount of dark 
energy density. In particular we consider a ACDM model character- 
ized by a large value of the dark energy density, JId e =0.9, which 
we refer to as L-ACDM (here L meaning large Qde value), and a 
cold-dark matter dominated cosmology, SCDM*, with Q,de ~ 
or equivalently Sim = 1 (the * symbol is to remind that the other 
model parameter values differ from the SCDM usually considered 
in the literature^). In Table|2]we quote the toy-model parameters, 
while all the other cosmological parameters (h, ris, erg, ..) are set to 

SCDM with Q-g = 0.51, h = 0.5 and shape parameter F = 0.5 as in 
Ijenkins et al] <1998l) 



Model noE a P(k) 



L-RPCDM 


0.74 


10 ACDM-W5 


L-ACDM 


0.9 


ACDM-W5 


SCDM* 





ACDM-W5 



Table 2. Toy-models characteristics; the values of the other parameters h, 
erg and Us are set to the ACDM-W5 values quoted in Table[T] 

the ACDM-W5 values. An important point is that given the same 
initial conditions at some early time (e.g. at recombination) these 
models predict different matter power spectra at z = 0. However, 
we would like to investigate deviations from a universal behaviour 
of the mass function independently of the present form of the mat- 
ter power spectrum. Therefore we artificially force these models to 
have the matter power spectrum at 2 = of the ACDM-W5 ref- 
erence cosmology. From a practical point of view this means that 
when running the N-body simulations, for each model we gener- 
ate the initial conditions such that the linearly extrapolated matter 
power spectrum at z — (obtained by using the growth factor 
(z) specific to each model) coincides with that of the ACDM- 
W5. Since D^{z) is model dependent, it implies that the various 
models will have the same initial power spectrum, but a different 
initial redshift (see Sect[3]and Sect|4]for technical details about this 
point). We want to stress that such toy-models are not intended to 
be compatible with observations, c ontrary to the "real istic models" 
considered in the previous paper ( lAlimietal.ll2010l) which were 
calibrated against CMB and SN la data. Again, our aim here is to 
perform a physical study of the cosmological dependence of the 
mass function. Nonetheless the conclusions drawn from this study 
will be extended to more realistic cosmological models as well in a 
forthcoming paper. 

In Fig.|2]we plot the Hubble rate H{a), the acceleration pa- 
rameter q{a) and the linear growth factor D'^ (a) /Dq for different 
cosmologies including our toy-models. The Hubble rate and the 
growth factor are normalized to that of the ACDM-W5 cosmology. 
In the upper panels we plot RP models with a = 0, 1, .., 10; the 
L-RPCDM model is shown as triple-dot-dashed line. We may no- 
tice that the Hubble rate increasingly deviates at higher redshifts 
from that of the reference cosmology for larger values of a (curves 
bottom to top in the upper left panel). Similarly, in such models 
the evolution of the acceleration parameter shows that the expan- 
sion is less decelerated during the matter dominated era and less 
accelerated during the dark energy dominated era (see upper cen- 
tral panel). Besides the acceleration starts later for larger values 
of a, such that above some values (corresponding to dark energy 
models with equation of state w > — 1/(30db)) the accelera- 
tion never takes place. As dark energy dominates earlier (for in- 
creasing values of a) the deviation of the linear growth rate with 
respect to the ACDM-W5 case is larger at high redshift, with dif- 
ferences up to a factor 2 at a = 0.1 (see curves bottom to top in 
the upper right panel). In the middle panels of Fig. |2] are shown 
ACDM models (equivalent to RP models with a — 0) with differ- 
ent amount of dark energy density. The various curves correspond 
to Qde = 0, 0.1, .., 1 (curves bottom to top). The L-ACDM and 
SCDM* models are plotted as short-dashed and dot-dashed lines 
respectively. We can see here a trend which is similar to that of 
varying a in the RP models. However when increasing Qde the 
deviations on the Hubble rate are much larger compared to the pre- 
vious case. Moreover, since the dark energy equation of state here 
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Figure 2. Scale factor dependence of the Hubble rate (left column panels), acceleration parameter (middle column panels) and linear growth factor (right 
column panels) for the RP models with different values of o = 0, 1, .., 10 (upper row panels), for ACDM with different values of VloE = 0, 0.1, .., 1 
(central row panels) and for the WMAP models (lower row panels). In all panels the solid line corresponds to the ACDM-W5; in the upper panels the L- 
RPCDM model corresponds to the triple-dot-dashed line, while in the central panels the L-ACDM and SCDM* models correspond to the short-dashed and 
dot-dashed lines respectively. In the bottom panels ACDM-Wl and ACDM-W3 models corresponds to the dotted and long-dashed lines respectively. 



is fixed to the cosmological constant value (q = or w = — 1), for 
decreasing values ofQ.DE the acceleration occurs later, and even- 
tually never takes place for sufficiently low values ofQ.DE (corre- 
sponding to matter dominated cosmologies). On the other hand we 
may notice that the linear growth factor deviates from that of the 
reference cosmology both at low and high redshifts (middle right 
panel). Finally in the lower panels are shown the ACDM- WMAP 
models. As we can see these have very similar behaviours both 
at the background and linear level. Henceforth our toy-models are 
characterized by a cosmic evolution which differs from that of stan- 
dard vanilla ACDM cosmology. Deviations on the Hubble rate are 
up to 40% for the L-RPCDM model, 50% for the L-ACDM and 
95% for SCDM* respectively. Also the redshift dependence of the 
linear growth factor differs from that of the reference cosmology 
with larger deviations at higher redshifts, up to a factor 2 at a = 0.1 
for the L-RPCDM model. In the case of L-ACDM and SCDM* 
such deviations are present both at high and low redshifts and up 



to 30% level respectively. If the universality hypothesis holds then 
such differences should not leave any imprint in the halo mass func- 
tion. 



2.3 Cosmological toy-models: spherical collapse 

The collapse of a spherical matter overdensity embedded in an 
expanding Friedmann-Lemattre-Robertson- Walker (FLRW) back- 
ground is the simplest model to describe the evolution of a density 
perturbation throughout the non-linear phase of gravitational col- 
lapse. For a given cosmology it allows us to estimate the time of 
collapse of an initial density perturbation, as well as the value of 
the overdensity at the time of virialization (i.e. assuming that the 
perturbation virializes at some time before collapse). Formally, the 
spherical collapse equation describes the evolution of a matter over- 
density 5m with a top-hat profile in a sphere of radius R as derived 
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Figure 3. Evolution of 5c (left panel) and A^ir as a function of the redshift of collapse for the different cosmological models, the line coding is the same as 
in Fig. [2] Note that 5c and Delta^i,, for L-RPCDM does tend towards 1.686 and 178 but at very high redshift. 



from the Raychaudury equation and which reads afl: 



R 
R 



4ttG 



(p + 3p), 



(13) 



where p and p are the total energy density and pressure within the 
sphere. Here we treat dark energy as a homogeneous component 
presents only at the background level. This is required for consis- 
tency with the fact that our N-body simulations do not explicitly 
account for the presence of dark energy perturbations, except for 
the linear regime of the realistic model simulations discussed in 
jAlimietal. 120101) . Hence the halos detected in the simulations are 
effectively dark matter overdensities embedded in a FLRW back- 
ground in the presence of a homogeneous dark energy. Since we 
aim to interpret the properties of these halos using the spherical 
collapse model, we solve the model equations in a similar cosmo- 
logical setup. 

In such a case the energy density and pressure in the spherical 
region are given by p = + poE and p — wdePde, where 

= pm(l + Sm) is the matter density in the sphere of radius 
R, while pm and poE are the background matter and dark energy 
density respectively, with w being the dark energy equation of state. 
The background densities evolve according to 



and 



Pm + 3-pm = 0, 

a 



pDE + 3(1 + w)-pnE = 0. 
a 



(14) 



(15) 



For the quintessence models considered here the dark energy den- 
sity and equation of state depends are given in terms of the scalar 
field potential and kinetic energy poE ~ 4>^/2 + V{(j>) and 
■"^ ~ \J?'^2\}f'}\ ' ths scalar field <j!> evolving according to the 
Klein-Gordon equation which we solve numerically. The evolution 



^ For a recent study of the spherical collapse model in dark energy cos- 
mologi es and a fully relativistic derivation of the Raychauduiy equation we 



of the local matter density inside the sphere is given by 

R 

Pm + 3;^Pm = 0. (16) 

Integrating Eq. ( I16t and substituting the definition of p%^ in terms 
of the background matter density, we find the relation between 5m 
and R at time t. 



i+.m^(i+c)(^)7|" 



(17) 



where 5^, a; and Ri are the matter overdensity, scale factor and 
radius at the initial time ti. Finally using Eq. |[17) we can rewrite 
Eq. l ll3t in terms of a seco nd order differential equation for the 
spherical matter overdensity jNunes & Mota ll2006l : lAbramo et al. I 
2007, see also): 



+ 2// 5 

T7 



47rG'pm5m(l + (5m) + 



4 5^ 



(18) 



refer to lCreminelU et al. I j2009l) 



3 1 + 5m ' 

Equation l |18t is a non-linear ODE, which at the linear order re- 
duces to the standard equation for spherical linear matter pertur- 
bations i n the presence of a hom ogeneous dark energy component 
(see e.g. ICreminelli et al. Il2009l in the inhomogeous case). Notice 
that dark energy affects the matter spherical collapse only because 
of its influence on the background dynamics through the friction 
term appearing in Eq. l |18t . 

Past works have solved the spherical collapse using Eq. l |13K 
however we prefer to work with Eq. il8t which we solve numeri- 
cally as an initial conditions problem rather than a boundary value 
one as in the case of Eq. l |13t . This allows us to relax some the 
assumptions used in the literature, such as estimating the time of 
collapse tc assuming that tc — 2tta, where tta is the time of turn- 
around, i.e. when _R = after an initial phase of radial expansion. 
Besides by solving the linearized form of Eq. l llSt for the same set 
of initial conditions that lead to collapse at tc (or equivalently Zc), 
we are able to directly infer the value of the linear density contrast 
at the time of collapse, 5c, without the need of using semi-analytical 
formula for the linear growth factor. In order to numerically solve 
Eq. l llSt , we first rewrite it in terms of the redshift variable z (i.e. 
It ~ "SzW' hsrsafter ' — d/dz) and set initial conditions with the 
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perturbation 5^ initially starting in the Hubble flow (Ri/Ri — Hi, 
i.e. Sl^ — 0). At the beginning the system evolves linearly, but 
as the quadratic terms in Eq. J18t become non-negligible the col- 
lapse enters in the non-linear evolution which ultimately leads to 
a diverging non-linear matter density contrast, 5™, — >■ oo, at the 
redshift of collapse Zc (corresponding to the radius of the spherical 
region R 0). Since our goal is to determine the linear density 
contrast for a perturbation collapsing at Zc = and more gen- 
erally at a given Zc, we have implemented a numerical algorithm 
which iteratively searches for the initial density perturbation value 
Sin for which 5m diverges at the input redshift value Zc. Formally 
5,n{zc) ~ OO, but numerically it is not possible to verify such a 
condition exactly, hence we use an additional algorithm, similar to 
that developed in [London & Flannerv 1 ( Il982l) , to determine for a 
given level of accuracy when Sm reaches the singularity. We fix the 
accuracy to Azc — 10~^, then the algorithm alts the search for 5^ 
as soon as a divergent solution for S,n is found within Azc from the 
specific Zc value of interest. At the same time we calculate on flight 
the turn-around condition S'^ ~ and determine the correspond- 
ing redshift zta- Furthermore we numerically solve the linearized 
form of Eq. dlSt using 5^ ~ and the 5^ value previously deter- 
mined for which collapse occurs at Zc to finally infer the value of 
the linear density contrast at the time of collapse. Sc. 

From the spherical collapse dynamics we can also calculate 
the redshift and the non-linear overdensity value at the time of 
virialization. In the spherical collapse model, this can be imple- 
mented only as an external condition by requiring that at the time 
of virialization the kinetic energy of the system T and the grav- 
itational potential energy U s atisfy the viria l condition T^ir = 
{R^ir/2){dU/dR)^ir (seee.g. |Peeblesll993h . Then using energy 
conservation, Tvi,- + U(zvir) ~ const., one can infer the redshift 
of virialization by solving the algebraic equation 

R DU 

Y^(z.«r) + (7(2:«zr) = C/(2ta), (19) 

where we have used the fact that at turn-ar ound the kinetic energy 
of the system vanishes, Tta ~ 0. We refer to lMaor & Lahav 
for an explicit form of the gravitational potential energy in terms of 
the matter overdensity and dark energy density for different cosmo- 
logical set up, including SCDM, ACDM and homogeneous dark 
energy models. Using the numerical computation previously de- 
scribed we numerically solve Eq. J19t to derive z^jir and compute 
the virial overdensity A„ir = Pm(2t>ir)(l + 5'^'')/ pm{zc)- How- 
ever we would like to remind the reader that the use of Eq. l |19b is 
rigorously justified only in the SCDM and ACDM models. This is 
because in the presence of a homogeneous dark energy component 
for which the background density evolves in time at a rate different 
from that of the background matter component, energy is not con- 
served within the spherical overdensity region. De spite several at- 
tempt s to account for such a loss of energy (see e.g. lMaor & Lahav] 
2005), we still lack a full relativistic calculatior0 (see discussion in 
Crem inelli et al. 2009). Thus similar to previous analysis the virial 
overdensity value which we infer for the L-RPCDM model should 
be considered only as approximative. Given this state-of-art com- 
putation we cannot do better. 

In Table[3]we list the values of Sc and Avir at Zc = for the 
different cosmological models. We recover the standard SCDM* 
values Sc = 1.686 and Avir = 178; the values corresponding 
to the various AC DM models are consistent with those found in 
lEke et al. I ( Il996h , similarly the values of the L-RPCDM model 

* We thank Jorge Norena and Filippo Vemizzi for pointing this to us. 



Model 


Sc 




ACDM-W5 


1.673 


368 


ACDM-W3 


1.672 


387 


ACDM-Wl 


1.674 


344 


L-RPCDM 


1.638 


436 


L-ACDM 


1.665 


708 


SCDM* 


1.686 


178 



Table 3. Values of the linearly extrapolated critical density threshold Sc 
and virial overdensity A^ir at = predicted by the spherical collapse 
model. Note that ACDM-WMAP cosmologies have very similar values 
for the collapse and virialization parameters. The L-RPCDM model dif- 
fers mainly in the value of Sc, while the effect of varying the dark energy 
density as described by the SCDM* and L-ACDM models primarily affects 
the virial overdensity value. 



are compatible with those quoted in lMainini et al. I (l2003h . As ex- 
pected, the ACDM-WMAP cosmologies have very similar collapse 
and virialization parameters, this is coherent with the fact that their 
expansion history and linear growth evolution are very similar as 
well. In contrast the L-RPCDM model predicts the lowest value, 
Sc — 1.638, which in the framework of the Press-Schechter for- 
malism means that in such a model structures form earlier. The 
spherical collapse prediction of Sc for L-ACDM and SCDM* gives 
Sc = 1.665 and 1.686 respectively, only a few percent lower and 
higher than those of ACDM-WMAP models. In contrast the former 
models predict very different values for the virial overdensity (708 
and 178 respectively). Since a larger virial overdensity implies a 
more compact object, we expect that virialized objects tend to be 
more compact as the amount of dark energy increases. 

In Fig. [3] we plot the evolution of Sc (left panel) and A^jir 
(right panel) as a function of the redshift of collapse. We can see 
that both the collapse threshold and the virial overdensity is closed 
to the SCDM values with large Zc- This is expected since at higher 
redshift the collapse occurs in a matter dominated universe. The 
redshift evolution of Sc in L-RPCDM shows the largest devia- 
tion with respect to the trend of the other models, with a much 
slower convergence toward SCDM. This is because dark energy 
starts dominating earlier in L-RPCDM than the other cosmological 
models. In contrast the evolution Avir(zc) has the largest deviation 
for the L-ACDM model, with Avir rapidly decreasing towards high 
redshifts due to the fact that in such a model the acceleration occurs 
earlier (see Fig.|2]l. 

The spherical collapse model remains a very simplistic de- 
scription of the gravitational collapse of structures in the universe. 
Typically these are not isolated, or spherical, and their velocity dis- 
persion and angular momentum are certainly not negligible. Nev- 
ertheless the model captures some fea t ures of the non-linear phase 
of collapse (see for instance FvalageasI ^2009 ) for analytical results 
in the rare events limit). It can therefore guide us toward a better 
understanding of the formation and evolution of dark matter halos 
as detected in N-body simulations. 
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3 N-BODY SIMULATIONS 
3.1 Simulation sets 

The numerical se t up is tiie same as ttie one presented in 
lAlimi et al. and we refer tiie interested reader to tiie pre- 

vious article for a detailed description of the numerical codes 
used in this series of papers. The N-body simulations are per- 
forme d using the RAMSES-AMR (Adap tative Mesh Refinement) 
code jTevssierll2002l :l Rasera & Tevssier I 12OO6.) based on a multi- 
grid Poisson solver. T he matter power sp ectra are computed us- 
ing the CAMB code (^Lewis et al. I l2000t) . then Gaussian initial 
conditions using the Zel'do v ich ap proximation are generated with 
MPGRAFIC ( Prunet et al. i ( l2008h '). For the toy-models we use 
the LCDM-W5 power spectrum (see previous sections). For each 
model the various cosmological variables (a{t),H{t), etc..) are 
pre-computed with an independent code and stored into tables, 
these are subsequently interpolated to input the various time depen- 
dent quantities into appropriately modified version of MPGRAFIC 
and RAMSES. We have used the same phase for the initial condi- 
tions of the various simulations, more specifically the white noise 
from the Horizon Project. We have performed a total of 15 simula- 
tions at the "Centre de Calcul Recherche et Technologie' ' (CCRtQ). 
Our longest simulation run for 350 hours of elapsed time on 64 
cores. In Table|4]we list the characteristics of the 15 cosmological 
simulations covering 6 cosmological models at various scales. 

Since we expect deviations from universality to manifest in the 
high-mass tail of the mass function, our box lengths of interest are 
648 h-^Mpc and 1296 h'^Mpc. In the case of the WMAP ACDM 
cosmologies we have also performed three simulations with a box 
length of 162h^^Mpc for comparison with the mass functions at 
intermediate mass ranges which have been estimated in previous 
studies. In Fig.|4l we illustrate the dynamics of the RAMSES AMR 
code with a zooming sequence of images of the dark matter parti- 
cle distribution in the ACDM-W5 at z — from cosmological to 
galaxy group scale (clockwise from top to bottom 648 h~^Mpc, 
162h"^Mpc, 40h'^Mpc and lOh^^Mpc). Because of the large 
simulated volume and the high spatial and mass resolutions, these 
simulations provide a robust statistics of well resolved halos. 



3.2 Halo finder and mass function 

The halos in the simulation boxes ar e identified using the Friend- 
of-Friend halo finder (FoF) jPavis et al. 1985). This algorithm de- 
tects halos as group of particles characterized by an intra-particle 
distance smaller than a given linking length parameter b. The al- 
gorithm runs over a list of particle coordinates in the box, firstly 
it regroups all those which are within distance b of an initial parti- 
cle. Then it moves to the next particle of the group and reiterates 
the same selection until no new neighbours are found. At this point 
the algorithm moves onto the next untagged particle repeating the 
above procedure to detect a new group (halo). 

Halos ca n also be detected usi ng the Spherical Overdensity al- 
gorithm (SO. lLacev&CoIelll994h . which provides halo mass es- 
timations similar to that performed observationally, moreover den- 
sity profiles are directly given by the algorithm. However the halo 
detection is restricted to spherical geometry, halos can overlap, and 
the last particles of the list could merge into some unphysical halos. 
In contrast under the FoF algorithm all particles belong to halos as 
in the case of the halo model, and it is applicable to non-spherical 



5 www-ccrt. cea.fr 



Model 


162h-iMpc 


648h-iMpc 


1296h-iMpc 


ACDM-W5 

Zi 

nip (h-^M©) 
Ax (h-^kpc) 

■^max 


93 

2.28 X 10^ 
2.47 
7 


56 

1.46 X 10" 
19.78 
6 


41 

1.17 X 10^2 
39.55 
6 


L-RPCDM 

nip (h-^M©) 
Ax (h-^kpc) 

■^max 


- 


136 
1.46 X 10" 
9.89 
7 


97 

1.17 X 10" 
39.55 
6 


L-ACDM 

Zi 

nip (h-iM0) 
Ax (h-ikpc) 

■^max 


- 


71 

5.61 X IQio 
9.89 
7 


52 

4.50 X 10" 
39.55 
6 


SCDM* 

Zi 

nip (h-iMg) 
Ax (h-lkpc) 

■^max 


- 


42 

5.61 X 10" 
19.78 

6 


30 

4.50 X 10^2 
79.10 

5 


ACDM-W3 

Zi 

nip (h-iM0) 
Ax (h-^kpc) 

■^max 


84 

2.11 X 10^ 
2.47 
7 


52 

1.35 X 10" 
79.10 
4 


38 

1.08 X 10^2 
158.20 
4 


ACDM-Wl 

Zi 

nip (h-iM©) 
Ax(h-ikpc) 

^max 


110 
2.55 X lO'' 
2.47 
7 


65 

1.63 X 10" 
19.78 
6 


47 

1.31 X 10^2 
39.55 
6 



Table 4. Parameters of the 15 N-body simulations for the various cosmolog- 
ical models: Zi is the initial redshift, mp is the mass of the particle, Ax the 
comoving resolution and I^sm. is the maximum refinement level. Each sim- 
ulation has 512'^ particles, with a 512'^ grid on the coai'se level and evolved 
down to 2 = 0. All the simulations share the same realization of the initial 
conditions (namely the Horizon Project white noise), and start at the same 
level of rms density fluctuation (~ 0.05) at the scale of the resolution of the 
coarse grid. Our refinement strategy consist in refining when the number of 
particles in one cell is greater than 8. 



particle groups as well, though it has tendency to overlink bridged- 
halos. Therefore FoF and SO are complementary algorithms. Here 
we use FoF since by relaxing the assumption on spherical symme- 
try this allows for a better physical study of the collapsed structures 
in the simulations (as we shall see hereafter). 

As can be seen in Fig. [4] halos are much more complicated 
than spherical isolated objects with a clearly defined radius. Halos 
in numerical simulations are non-spherical, interacting with their 
environment, contain lots of subhalos and, moreover the density 
can vary rapidly. This makes the definition of their boundary some- 
what arbitrary. The SO halo definition is based on the enclosed 
overdensity A in a given spherical region. This implies that it is 
rather straightforward to run a SO halo finder with threshold given 
by Auir from the spherical collapse. Similarly, it is possible to run 
a FoF halo finder with a linking length parameter h = b^ir, where 
bvir corresponds to the value of a spherical overdensity A„i,,, pro- 
vided that a relation between the two is assumed. A practical con- 
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4 NUMERICAL CONSISTENCY CHECKS 

In this Section we present a detailed study of various source of 
numerical systematic effects. As result of this analysis we are able 
to control numerical errors within few percent level, thus validating 
our final results at the same percentage level. 




Figure 4. 2-dimensional projection of the dai'k matter density at z=0 from 
the ACDM-W5 simulation in the 648 h^^Mpc box. The four images rep- 
resent a zooming sequence (clockwise from top to bottom) of a factor four 
for each image from the full box size (upper left) to a 10 h^^Mpc scale 
(bottom left). 



versio n formula is given by A/178 ^ (0.2/6)3 jCole & Lacev I 
(I1996')). nevertheless such a conversion is very approximative. A 
thorough comparison between the FoF and S O algorithms and the 
accuracy of conversions ha s been performed in lLukic et al. I ^20091) 
(see also lAudit et al. |[l998h . Furthermore since bvir depends on an 
approximate treatment of the non-linear collapse, its use as a link- 
ing length parameter, introduces built-in assumptions on the mass 
function obtained with FoF(bvir), which may hide the non-linear 
contributions associated with the halo formation and virialization 
process. For these reasons we start by considering the FoF halo 
finder with a constant b, which we set to 6 = 0.2. As we will show, 
our conclusion about universality will be independent of the exact 
value of b, remaining valid in the range 0.1 to 0.3 (corresponding 
to Avir varying roughly from 1424 to 53). 

One last point concerns the numerical definition of the halo 
mass function, which reads as 



TV 



dn{M) _ J 

dlog(M) ^ U Alog(M) '' 



(20) 



where N is the number of halos in a logarithmic mass bin between 
log(M) - Alog(A/)/2 and log(A/) + Alog(M) /2, with Alog(M) 
the bin width and L the comoving box length. Our default choice 
of the bin width is very conservative, AM/M ~ 0.2. Moreover 
when measuring the mass function we focus on mass ranges which 
corresponds to halos containing at least 350 particles and Poisson 
shot noise below the 10% level. As we will discuss in the next 
Section, such conservative assumptions are necessary to limit the 
effect of mass resolution, bin size and Poisson noise, which can be 
important source of systematic errors in the high-mass tail. 



4.1 Initial conditions 

We have performed a series of tests on the initial conditions 
generated with MPGRAFIC. Firstly, we have checked that the 
power spectrum of the initial particle distribution is consistent 
within a few thousandth with the input spectrum from CAMB 
near the Nyquist frequency (at this high frequencies the Poisson 
noise is minimal). In addition, we have compared the MPGRAFIC 
power spectrum with that extracted from the GRAFIC code for 
the same white noise. After turning off the Manning filter subrou- 
tine in GRAFIC we have found good agreement even at low k 
( iBertschinger II 199511200 ih . The Hanning filter suppress the high k 
modes, hence potentially causing unphysical numerical artifacts on 
the generated particle distribution even in the linear regime. There- 
fore previous studies that have inadvertently used the Hanning filter 
may have underestimated the mass function, which is known to be 
sensitive to the gravitational dynamics on all scales. In our case, the 
use of the Hanning filter, led to a roughly 10% suppression of the 
mass function over the full range of masses tested in our simula- 
tions (10^* — lO^^h^^Mo). Therefore, we preferred not to use this 
filter. 

A further test concerns the choice of the initial redshift 
of the simulations. This can be responsible for spurious ef- 
fects due to transients from initial conditions when using the 
Zel'd o vich approximation (Reedetal. 2003, 2007;^ Crocce et al. 



200d iLukic et al. 1 12007|). Several works iCrocce et al. I fcoO^ . 
Tinker et al. I ( I2OO8I) and lcrocce et al] ilOm suspect that a very 
low initial redshift combined with the Zel'dovich approximation 
might be responsible for the discrepanci es between the i r estim ated 
mass functions and those inferred by llenkins et al. ] i200lh and 
IWarren et al. To be as conservative as possible we have 

chosen higher initial redshifts Zj compared to previous s tudies 
dSheth et al. II2OOII ; llenkins et al. Il200ll : lwarren et al.1l2006h . This 
is realized by imposing that the standard deviation of the density 
fluctuations smoothed at the scale of the coarse grid is 
^^^coarse-j ^ q gg ^-jj^ jj^j^ ^-jjoi^-g for the different cos- 

mologies and simulations boxes vary in the range 30 — 110, which 
ar e much higher then f or instance the initial redshifts considered 
m IWarrenetal7l l l2006h which range from 24 to 34 depending on 
the box lengths, and for which lCrocce et al~l ( 12006) have estimated 
the suppression of the high-mass tail of the mass function to be of 
order 10 - 15%. 



4.2 N-body code 

We have run a series of tests to check the accuracy of the RAMSES 
code as well as the modifications that have been implemented to 
account for the quintessence scenario. Firstly, we have verified that 
for the various models the ratio of the measured power spectrum to 
the linear prediction from CAMB is constant at the percent level on 
the large linea r scales (k::::i0.01-0 . 1 h/Mpc) of the simulations at all 
redshifts (see lAlimi et al. ir201Cll) . Furthermore, we have checked 
that by reducing the integration time steps the results remain un- 
changed, thus validating our implementation of quintessence mod- 
els in RAMSES. Secondly, we have compared the mass function 
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from the RAMSES simulation of tlie ACDM-W3 cosmology withi 
box lengtii 648 h~^Mpc and 256"^ particles against a simulation 
with identical characteristics and same see d for the initial con- 
ditions obtaine d using the GADGET code jSpringel et al. lIzOOll : 
ISpringellllOOSl) . We find that the two mass functions are in very 
good agreement (few percent level) over the range of masses corre- 
sponding to halos with m ore than 350 particle s and less than 10% 
Poisson noise. Using the Warr en et al. I ( 1200^ correction does not 
improve noticeably the inferred mass fun ction below 350 parti cles. 
In principle the halo mass correction bv I Warren etdl j2006h de- 
pe nd on code and/o r cosmology. Another point concerns the results 
of iHeitmann et al.1 (2005, 2008), the authors have computed the 
mass function with various codes using the same set of initial con- 
ditions. They found a s catter of order 10% in the high mass range, 
the HOT code used in I Warren et all bOOd) seems to give lower 
mass functions than other codes such as GADGET by a factor of 
about ten percent at ~ 10^^ h'^M©. 

4.3 Particles, box length, mass and force resolution 

Another set of tests concerns the sensitivity of the results to var- 
ious simulation characteristics. First we have checked the stabil- 
ity of the results when varying the number of particles from 256'^ 
up to 1024'^, while keeping the same phase of the initial condi- 
tionfl Assuming 350 particles as a minimum for the halo mass, 
we find that the mass functions are consistent to better than 5% at 
all masses. Secondly, we have evaluated the influence of varying 
the maximum level of refinement in RAMSES from 6 to 4. This 
parameter changes the spatial resolution of the simulations, how- 
ever we found no effect on the mass functions. Besides in our sim- 
ulations the refinement level evolves freely, and never reaches the 
maximum allowed level. 

Another important consistency check concerns the coherence 
of the mass functions as measured from simulations with different 
box lengths. Discrepancies may indicate effects due to either mass 
resolution or finite volume. With our conservative choice of a min- 
imum of 350 particles per halo and 10% level of Poisson noise, 
the coherence in the mass range of interest is better than 5%, as it 
can be seen in the left panel of Fig.|5] where we plot the difference 
between the mass function measured in ACDM-W5 simulations at 
z = with 512'^ particles and box lengths of 648 h~^Mpc and 
1296 h~^Mpc respectively, and the Sheth-Tormen parametriza- 
tion, Eq. ((S)), calibrated over the whole set of ACDM-W5 simu- 
lations. The vertical dotted lines identify the mass intervals where 
the Poisson noise is below 10% and halos have at least 350 parti- 
cles (the right interval corresponding to 1296 h~^Mpc simulation 
box and the left one to 648 h^^Mpc). We can see that the scatter 
of the points around the zero residual is well within the 5% level. 
We did not investigate finite volume effects, or cosmic variance 
specifically. In fact in our largest box (1296 h~^Mpc), the cos- 
mic variance error on the largest halo mass es is likely to be dom- 
inated by Poisson shot noise as shown in jHu & Kravtsov ll2003l : 
ICrocce et alj|201ol) , hence it should be within the 10% level set 
by our requirement for the Poisson noise. We have also checked 
the effect of changing the phase of the initial conditions. In the 
right panel of Fig. |5] we plot the residuals of the measured mass 
function in ACDM-W5 simulations at z = with box lengths of 

^ The results of these simulations with I billion particles are pai1 of the 
Dark Energy Universe Simulation Series (DEUSS) and will be presented in 
an upcoming paper. 



1296 h~^Mpc and 512^ particles with two different phases, with 
respect to the ACDM-W5 best fitting formula. In the mass range of 
interests the differences between the residuals are well within the 
Poisson errorbars. 



4.4 FoF code & Mass Binning 

We have compared the results fr om our halo finder code to those 
obtained using "Halomaker" FoF jTweed et al. |2009|) and found no 
difference in the measured mass function. We have also performed 
a study of the influence of varying the bin-width and the binning 
strategy on the mass function. Th ese tests are particul arly impor- 
tant, since some authors (see e.g. Ijenkins et al. II2OO1I) have used 
analytical corrections to remove effects caused by large mass-bins 
and smoothing functions. 

In Fig. |6] we plot the residual mass function of ACDM- 
W5 simulations a.i z — Q with 512^ particles and box lengths 
of 648 h~^Mpc and 1296 h~^Mpc respectively with respect to 
the Sheth-Tormen best fitting formula, for different bin-widths: 
AM/M = 0.1 (left panels), AM/M = 0.25 (central panels) and 
AM I'M = 1 (right panels). In the top panels we plot the case with 
averaged-mass bins, i.e. the mass of each point is the average mass 
of the halos in the bin, while in the bottom panels we plot the case 
with centered-mass bins, i.e. the mass of each point is the central 
value of the logarithmic mass bin. We can see that for both binning 
strategies increasing the bin-width reduces the Poisson noise and 
smooths the curves. However increasing the bin size is particularly 
delicate in the high mass tail especially for the averaged-mass bin- 
ning. In fact we can see that a larger binning {AM/M = 1) tends 
to lower the mass function intrinsic scatter. This is because as the 
bin-size increases, each bin has a larger number of halos thus lead- 
ing to smaller Poisson errors on the other hand in the high-mass tail 
the mass function drops steeply and the use of averaged-mass bins 
tends to lower its value. As noticeable from Fig. (6] such an effect 
is much weaker in the case of centered-mass bins which we adopt 
hereafter, and for the study of the universality of the mass func- 
tion we consider centered-mass bins with width AM/M — 0.2, 
a compromise between t he Poisson noi s e and the number of bins. 
It is worth noticing that Ijenkins et al. I bOOlh use large bins and 
a gaussian smoothing (of rms 0.08 in In (M) corresponding to a 
width of 0.37 in our units) which is supposed to strongly raise the 
curve. An analytical correcting factor is then applied to recover a 
proper estimate of the mass function. We find this procedure to be 
quite uncertain. Again as seen in Fig [6] using large bins tends to 
lower the mass function unlike expectations, and such an analytical 
correction might cause an underestimate in the high mass end at 
the 10% level. In such a case, it is safer to use smaller bins such 
that deviations remain within 5% (as in our case) and not rely on 
analytical corrections. In the light of these various tests, our preci- 
sion on the absolute FoF mass function is expected to be of order 
5 — 10%. However, we want to stress that the precision on the rela- 
tive mass function is even higher (typically less than 5%) since for 
all our simulations we use the same phase for the initial conditions, 
and thus systematic errors cancel out. The Poisson error bars are 
consequently not relevant for quantifying the error on the relative 
mass functions, since the scatter of the mass functions in two dif- 
ferent cosmologies is strongly correlated. Moreover from the above 
analysis we can see that intrinsic fluctuations of the measured mass 
functions for a given cosmology are much smaller than the Poisson 
error bars. 
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Figure 5. Left: residual between the measured mass function and our fitting formula for two ACDM-W5 simulations with box lengths 162 h^^Mpc and 
648 h^^Mpc at z=0. The number of particles is 512'^, the minimal number of particles per halo is 350 and the Poisson noise is less than 10% (vertical lines). 
This illustrates the coherence between the simulations. Right: As in the previous figure for two simulations with different realisations of the initial conditions. 
The box length is 1296 h^^Mpc and the number of particles 512^. 
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Figure 6. Effect on the mass function of varying the bin width and binning strategy. The plots show the residuals between the measured mass functions 
and our fitting formula Eq. i22\ as a function of ln{<T~^) for two ACDM-W5 simulations with box lengths 648 h~^Mpc and 1296 h~^Mpc at z=0. The 
minimal number of particles per halo is 350 and the Poisson noise is less than 10%. Left: bins width AM/AI = 0.1. Middle: 0.25. Right: 1. Top raw: 
"average-mass" bins. Bottom raw: "Centered-mass" bins. Using too large bins decrease the mass function, especially for the "average-mass" bins. To be as 
conservative as possible we use "centered-mass" bins of width AM/M ~ 0.2. 
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Figure 7. Mass functions dn{M)/dlog{M) for the ACDM-WMAP cos- 
mologies at z=0 (upper curves), and z=l (lower curves) from 162 h~^M0, 
648 h~^MQ and 1296 h~'^MQ box lengths simulations. Halos are de- 
tected with FoF(fe = 0.2). The dotted line is our best fit ACDM-W5 
Eq. (22). Crosses correspond to ACDM-Wl, stars for ACDM-W3 and tri- 
angles to ACDM-W5. The mass functions cover a range from Milky- Way 
size halos to cluster of galaxies. 
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Figure 8. Mass functions for the ACDM-WMAP cosmologies at z=0, in 
the / — ln((T^^) plane (as defined in the text) from the 162 h^^M©, 
648 h~^MQ and 1296 h~^MQ box length simulations. Halos are detected 
with FoF(f) = 0.2). The continuous line is the ACDM-W5 fit described in 
this section (Eq. l22t and the dotted line is the Jenkins et al. (2001) fit. The 
conventions are the same than in Fig. [T] Deviations from universality are 
below numerical precision. 



5 RESULTS 

5.1 WMAP-ACDM cosmologies and universality at z = 

In Fig. |7] we plot the mass function of the ACDM-WMAP cos- 
mologies at 2 = and 1 measured in 162h~^Mpc, 648h~^Mpc 
and 1296 h~^Mpc simulation boxes respectively. Notice the large 
mass range covered by the simulations, which spans from MiUc- 
Way to cluster size halos (10^^ - lO^^h'^Mo). In order to test the 
universality of the mass function it is preferable to work with the 
function 



M dn 



(21) 



where dn/dln{a~^) is the comoving number of halos per unit of 
natural logarithm of a~^. As mentioned in Section lZT] if the mass 
function is universal the selection function accounting for the non- 
linear collapse in the definition of / is independent of cosmology 
(i.e. independent of the density threshold). In such a case the func- 
tion / for different cosmological simulations should be identical to 
numerical precision. 

In Fig. [8] we plot the function / for the three ACDM-WMAP 
models at z = in the range —0.8 < ln(cr^^) < 0.7. We can 
see that the data points nicely overlap to numerical precision, this 
is quite remarkable given the fact that these models have differ- 
ent initial conditions with very different erg values. This implies 
that changing the initial conditions does not break the universality 
of the mass function in the ACDM cosmolo g ies, which is in agree- 
ment with the results of lJenkins et al. I j200lh : l^^en et al. I j2006l) . 
Nonetheless we find that the functional form of / differs from 
standard fitting formula inferred in previous an alysis jSheth et al. I 
1200 ll : I Jenkins et al. 11200 ll IWarren et al. II2006I) . As an example in 
Fig.[8]we plot a gainst the simulation data-points the fitting formula 
Eq. ^ given in jjenkins et al. l200lh . Large deviations with respect 



to this fit occurs in the high-mass end, instead we find a better fit 
by using the functional form of the ST formula, with different pa- 
rameter values calibrated on the ACDM-W5 model simulations. In 
particular by fixing 5c to the spherical collapse model prediction of 
the ACDM-W5 cosmology, Sc = 1.673, our best fitting function 
(using b=0.2) reads as: 



2a 



1/2 



1 + 



5c 



-2p 



(22) 



with A = 0.348, a = 0.695 and p = 0.1. 

We plot in Fig. |9] the residual of the measured function / 
against Eq. i22\ . we can see that the ACDM-WMAP cosmolo- 
gies have a universal mass function to 5 — 10% accuracy level. 
For comparison we also plot the original ST formula Eq. ([8} with 
5c = 1.686 (dashed line), and the Jenkins et al. one Eq. ^ (dot- 
ted line). The former badly reproduce the measured mass functions 
at ln((7^^) > —0.5, while the latter provide a good description 
only in the lower mass range (ln((T~^) < 0.4). In contrast at higher 
masses there are deviations > 20%, thus confirm ing previous re- 
sults by [Tinker et aTl ( l2008h : ICrocce et al.l l2oTdh . 

The fact that in ACDM-WMAP cosmologies the mass func- 
tion at z = shows a universal behaviour is not surprising, since 
from the considerations of Section[2jT]these models despite having 
different erg values and small differences in the other parameters, 
they have nearly identical expansion histories, linear growth evo- 
lution and spherical collapse model parameters. The spherical col- 
lapse model is only approximative we may expect small deviations 
from universality to be present also in these models, however it is 
likely that such deviations are within the numerical errors of our 
simulations. 

The universality of the mass function for cosmologies charac- 
terized by the same expansion history but different ag values (and 
slightly different values of the other cosmological parameters) im- 
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Figure 9. The residual between the fitting formula Eq. \22\ and the mea- 
sured FoF(fe = 0.2) mass functions at z=0 for the "WMAP ACDM cos- 
mologies". The conventions are the same than in Fig. [T] Halos come from 
162 h-^MfT) 648h-^Mf7) a nd 12 96 h-^MQ box length s simulations. For 
reference the Ijenkins et al 1 <200lh (dotted line) and the ISheth & Tormeni 
<1999l) : IShethetal.U200ll) (here in dashed line) fitting formula are also 
plotted. We recover the universality up to numerical precision (5 — 10%). 



plies that the position of the mass function in the plane / — ln(cr^ ) 
is independent of erg (within the accuracy of our simulations). This 
is an important point to keep in mind, as we will show in the next 
Section this will allow us to isolate cosmological dependent effects 
on the halo mass function, which are related to the record of the past 
expansion history, and also extend our conclusions on the limits of 
universality to models with different as values. 



5.2 Dark energy models and the limit of universality 

We now focus on the toy-model simulations. In Fig. [TO]we show 
the projected density maps from the 648 h~^Mpc simulation box 
at 2; = zoomed on 40h~^Mpc (left panels) and 1 Oh ~^Mpc (right 
panels) scales respectively. Since the simulations have the same ini- 
tial phase, in order to facilitate a visual comparison between differ- 
ent models we plot their density distribution in the same image with 
a different color coding for each model. In the top panels are shown 
the density maps for the ACDM-W5 (red) and L-RPCDM (green), 
while in the lower panel in addition to the ACDM-W5, we also 
show the L-ACDM (green) and SCDM* (blue). The differences in 
the top panels are indicative of the effects related to having a time 
evolving dark energy with 1^(2:) > — 1, while those in the bottom 
panels are associated to varying the amount of dark energy density. 
We can see that in both cases there are no apparent differences in 
the particle distribution on the 40h~^Mpc scale. In contrast differ- 
ences are clearly manifest on the 10h~^Mpc scale, where the halo 
concentration seems to differ from one model to another, and the 
outer parts of halos show different morphologies as well. In par- 
ticular, we may notice that in the L-RPCDM case, since f2,„ and 
P{k) are the same of the reference cosmology, differences in the 
dark matter distribution are unique signature of the non-linear pro- 
cess of structure formation. In such a case, it is hard to believe that 
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Figure 11. The residuals between the measured mass functions of L- 
RPCDM (diamonds) and the measured mass functions for ACDM-W5 (tri- 
angles) cosmology at z=0. Deviations from universality are clearly above 
numerical errors (< 5%) and correlated with linear growth history. Here 
the amplitude of the residuals is a signature of the different dark energy 
equation of state evolution in the L-RPCDM model. 
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Figure 12. The residual between the measured mass functions of L-ACDM 
(crosses), SCDM* (squares) and that of ACDM-W5 (triangles) cosmol- 
ogy at z=0. Deviations from universality are clearly above numerical errors 
(< 5%) and correlated with linear growth history of the models. Here the 
amplitude of the deviation is related to the different amount of dark energy 
of the models. 



the mass function should remain unaffected, that is to say universal 
in dark energy cosmologies. 

For a more quantitative comparison, we plot in Fig. [TT] the 
residual of the function /(cr) measured in L-RPCDM with respect 
to that of the ACDM-W5. As previously mentioned we only plot 
the mass range in which halos contain at least 350 particles and 
where the Poisson noise is within 10%, thus allowing us to con- 
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Figure 10. Projected overdensity maps for the ACDM-W5 and L-RPCDM models (top) and ACDM-W5, L-ACDM and SCDM* (bottom) on scale 
40 h~^A/pc (left) and 10 h~^Mpc (right) respectively at z=0. The views are centered onto two halos and extracted from simulations of box length 
648 h-^Mpc with the same realisation of the initial conditions and the same trg. Each cosmology is represented with different color coding (but with 
the same intensity level). Top panels: ACDM (red) and L-RPCDM (green). Here the distinct green or red regions at small scales are non-linear imprints related 
to the different equation of state evolution between the two models. Bottom panels: ACDM (red), L-ACDM (green) and SCDM* (blue). Here the distinct 
colored regions on small scales are non-linear signature due to the different amount of dark energy of the models. 



trol numerical errors to better than about 5%. From this plot we 
can clearly see deviations from universality at the 10% level, hence 
well above numerical errors. Similarly in Fig. [12] we plot the case 
of L-ACDM and SCDM* . The former lies about 5% above the 
ACDM-W5 model, while the latter is roughly 10% below, again 
these are evidence of departure from universality. However such de- 
viations are not random, they are correlated with the linear growth 
history of the corresponding model. In fact comparing the evolu- 



tion of the growth factor of each model (see Fig. O with the cor- 
responding amplitude of the mass function residual, we may no- 
tice that the greater the deviation in the growth rate history and the 
larger is the deviation from universality at z=0. This is not the case 
for the ACDM-WMAP cosmologies which share nearly the same 
linear growth rate history and therefore their mass functions are 
universal to numerical precision. The physical origin of such devi- 
ations is quite intuitive. The cosmic structure formation is more 
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efficient in tiie past when tlie linear growth rate is greater. This 
is indeed the case for our toy-models. Even though by construc- 
tion they share the same amount of clustering at the linear level 
at 2 = 0, lots of dense halos were formed earlier and survive at 
lower redshift since they are decoupled from the cosmic expansion 
on the larger scales. The same effects are responsible for the im- 
prints wh ich h a ve bee n shown to affect the non-linear matter power 
spectrum Ma (2007); Alimi et al. ( 2010) and halo concentration 
iDolag et al. (2004); Wechsler et al. ( 2002) . To our knowledge, this 
is the first time that physical effects leading to a non-universal be- 
haviour of the mass function have been unambiguously shown. 



6 NON-UNIVERSALITY OF THE MASS FUNCTION 
AND THE NON-LINEAR COLLAPSE 

In this section we will focus on the non-linear processes which 
are responsible for the departure from a universal behaviour of the 
mass function. As discussed in Section IZTl deviations from uni- 
versality occur if and only if the non-linear collapse of dark mat- 
ter is cosmology and redshift dependent (or in other words if the 
relation between the linear and no n-linear growth of stru ctures is 
cosmology and redshift dependent Francis et al. I (|2009b')). In the 
Press-Schechter framework this is parametrized by the dependence 
of the mass function on the threshold density 5s- This suggests that 
we may gain a better insight into the origin of the non-universal be- 
haviour by explicitly accounting for the non-linear collapse model 
in the analysis of the measured mass functions. Here we use the 
spherical collapse model which has been described in Section |23] 
Furthermore we will plot residuals with respect to the mass func- 
tion measured from the SCDM* simulations, rather than that of 
the reference ACDM-W 5 cosmology. In fact theor e tical arguments 
l Efstathiouet"an ( Il988h ; Isianchard et al. I ( Il992l) ; iLacev & Cole I 
1 19941) ) suggests that for scale-free initial conditions (with a power 
law of slope — 1 < n < 1) in the Einstein-de-Sitter cosmology the 
structure formation should be universal (or self-similar) since no 
length or time scales (other than the one from the non-linear col- 
lapse itself) are involved. Though our power spectrum for the initial 
conditions is not a power law, the variations of the slope as a func- 
tion of k are very mild such that deviations from universality in the 
SCDM case should still be small. Moreover the spherical collapse 
model predict constant Sc and A^i,. values, and all our cosmologies 
converge to a SCDM-like behaviour at high redshift since the influ- 
ence of dark energy becomes negligible. Therefore residuals with 
respect to SCDM* should be more sensitive to the cosmological 
and redshift dependent signature of the non-linear collapse on the 
halo mass function. The mass function for the reference SCDM* 
used here for comparison is measured from the simulations of box 
lengths 648 Mpc/h and 1296 Mpc/h. As for the ACDM-W5 model 
we use a fitting formula with the ST form. We find for the SCDM* 
model the following parameter values: A — 0.350, 2 — 0.720 and 
p = 0.1 and 5c = 1.686. 

6.1 The role of the threshold of collapse 

Here we consider the effect of the critical density threshold (or 
equivalently the time of collapse) on the halo mass function as pre- 
dicted by the spherical collapse model. In the left panel of Fig. 1131 
we plot the residual of the function / for the L-RPCDM model 
with respect to the SCDM* case as a function of ln(CT~^). The am- 
plitude of the residual is of order 20% over the entire mass range. 
In the right panel of Fig.[T3]we plot the residual as a function of 
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Table 5. Standard deviations from a universal behaviour of the mass func- 
tion with and without accounting for the time of collapse as encoded by 5c 
predicted in the spherical collapse model. For comparison we also quote the 
values of 5c and A„ir for each model as in Table[3] 

ln((5»/(T), where 5* — Sc/1.686 is the critical density threshold of 
the L-RPCDM model rescaled to the SCDM* value. As we can see 
the deviation with respect to the SCDM* model is strongly reduced, 
more in the high mass end than at low masses. A similar trend oc- 
curs if we consider the residual of the L-ACDM and ACDM-W5 
models, which we plot in Fig. [14] as a function of \n{a^^) (left 
panel) and \n{5,/a) (right panel), with 5* the value of the criti- 
cal density threshold of each model rescaled to the SCDM* value. 
Again accounting for the collapsing time of halos as encoded in 5c 
reduces the residuals. For a more quantitative comparison we quote 
in Table[5]the standard deviation^ of the residuals of each cosmol- 
ogy as a function of a^^ and 5, /a respectively. We also quote the 
values of 5c and A„i,. at z = already discussed in Section [23] 
We may notice that deviations with respect to the SCDM* mass 
function are correlated with the difference in the value of 5c spe- 
cific to each model. In particular as 5c decreases with respect to 
the SCDM* value, the standard deviation increases up to ~ 20% 
for the L-RPCDM. Indeed accounting for the spherical collapse 
threshold reduces discrepancies between the halo mass functions 
by a factor ~ 1.4 for ACDM-W5 and L-ACDM, and a factor ~ 2 
in L-RPCDM. This clearly show that the density threshold plays a 
role in the cosmic structure formation and a general prescription for 
the mass function should include its dependence upon cosmology. 

From Table [5] we can also notice that the remaining residuals 
are still relevant, above numerical uncertainties, and not correlated 
with 5c. Moreover if we consider the L-RPCDM model, accounting 
for 5c improves only the high mass range and the deviations on the 
lower mass range can be hardly be attributed to 5c. This indicates 
that there is at least one more process responsible for the departure 
from universality, which will discuss in the next paragraph. 

6.2 Halo viriaUzation and redshift evolution of the mass 
function 

Confronting the level of improvement of the mass function residual 
as a function of 5t/a and the difference of the A^ir value of each 
cosmological model with respect to SCDM* quoted in Table|5]we 
notice a strong correlation between the two. The larger the differ- 
ence in the value of A^ir the greater the deviations of the function 
f{5t/o). Such a correlation is indicative of the fact that the viri- 
alization of halos does indeed play a role in determining the mass 
function and since the characteristics of this process are dependent 



The standard deviation used to quantify the differences with respect to 
the SCDM* mass function is defined as JY.i ^ log fi^/{N - 1) 
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Figure 13. Residual between the measured mass functions in L-RPCDM (diamonds) and SCDM* (squares) models respectively at z=0. Left: Residual for 
the mass functions is shown in the / — \n{l/a) plane. Right: Residual for the mass function in the / — ln((5«/o') plane, where (5* = 5c/1.686. Taking into 
account the different halo formation times through 5c strongly decreases deviations from universality in the high-mass end. 
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Figure 14. Residual between the measured L-ACDM (crosses), ACDM-W5 (triangles) mass functions and that of SCDM* (squares) at z=0. Left: Residual of 
the mass functions in the / — ln(l/cr) plane. Right: Residual of the mass functions in the / — ln(5* /a) plane, where 5* = i5c/1.686. Taking into account 
the different halo formation times through &c decreases deviations from universality in the high-mass end. 



upon cosmology and redshift, so are the deviations from universal- 
ity. Qualitatively this can be understood as follows. Let us consider 
two perturbations, one in SCDM* and another in L-ACDM, both 
collapsing and virializing at the same moments, then the spherical 
collapse model suggests that the resulting halos should be much 
more compact or "concentrated" in L-ACDM than in SCDM*. Be- 
cause we measure mass functions with a constant linking length pa- 
rameter (corresponding to a fixed overdensity), this would explain 
why the mass function is greater in L-ACDM than in SCDM* . 

Now, if we extend this argument to the redshift evolution of 
the virial overdensity, since for every cosmology Auir(2:) con- 
verges toward the SCDM* value at higher redshifts, then the mass 
function residuals must decrease with redshift in correlation with 



the behaviour of A„ir(z) specific to each model. In Fig.[T5] we plot 
the mass function residuals from the 648 h~^Mpc simulation boxes 
at 2 = 0, 0.25 and 1 for the SCDM* (upper left panel), ACDM-W5 
(upper right panel), L-RPCDM (lower left panel) and L-ACDM 
(lower right panelfl We can see that the deviations correlate well 
with A„ir {z) of each model. For the SCDM* the residual is consis- 
tent with a zero value, this is a non-trivial result which shows that 
clustering in SCDM models occurs in a universal manner as a func- 



° We hmit this analysis to the 648 h~^ box length and z < 1, since higher 
redshifts or larger box lengths would result in too few halos and larger Pois- 
son error bars in the reference SCDM* model. 
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Figure 15. Redshift evolution of the deviations from a universal behaviour. We plot residuals between the measured mass functions for SCDM* (squares, 
upper left panel), ACDM-W5 (triangles, upper right panel), L-RPCDM (diamonds, lower left panel), L-ACDM (crosses, lower right panel) with respect to 
the SCDM* mass function at z=0 (squares on the central line in the upper left panel). In each plot the left curve corresponds to 2 = 0, the middle curve to 
z = 0.25 and the right one to 2 = 1. The box length of the simulations is 648 h~^Mpc. As expected, all mass functions tend towards SCDM* at high 
redshift. 



tion of redshift. As shown in Section [23] the evolution of the virial 
overdensity in ACDM-W5 rapidly converges toward the SCDM* 
value and again in Fig.[T5]we can see that the mass function resid- 
ual rapidly vanishes as a function of redshift. However this also 
shows that the mass function as measured with FoF(6 = 0.2) for 
the ACDM-W5 is not universal for z < 1, where dark energy 
starts dominating th e cosmic expansion , which is in agr eement with 
recent findings bv iTinker et al. I j2008l) ; ICrocce et alj JioiO) . The 
models with the largest deviations of A„ir(2) are the L-RPCDM 
which converges towards the SCDM* value very slowly at high 
redshifts and L-ACDM which converges more quickly. In the lower 
panels of Fig. \\5\ we can see that this is indeed the case for the 
mass function residuals. Although errorbars become very large, we 
also checked that at higher redshift (for instance at z=2.3), the mass 



functions are all compatible with a null deviation at the same level 
independently of the cosmological model. 

These results clearly demonstrate that the virialization process 
also contribute to shaping the halo mass function in a cosmologi- 
cal and redshift dependent way. At this point we may ask ourselves 
whether accounting for the virialization in the measurement of the 
mass function may reduce the deviations from universality, sim- 
ilarly to the collapse threshold. After all we have detected halos 
using a constant linking length parameter, b — 0.2, and using in- 
stead a value of b as predicted by the spherical collapse may fur- 
ther reduce discrepancies of the mass function residuals. Hence 
we have run the FoF algorithm on the simulations at 2 = with 

S iarameter b^jr given by the conversion A„i,,/178 « (0.2/fe„ir)"^ 
Cole & Lacev 1 19961) , with A„i,. the value predicted by the spher- 
ical collapse model specific to each cosmology. As illustrated in 
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Figure 16. Residuals between the measured mass functions for the consid- 
ered cosmology at z=0 using FoF with b = byir (specific for each cosmol- 
ogy) from the spherical collapse and the mass functions for the SCDM* cos- 
mology. Conventions are the same as in previous graphs. Accounting for the 
virialization overdensity from the spherical collapse models strongly over- 
corrects the deviations from universality. Such deviations are larger than 
those one find s using a constant b for detection (which is consistent with 
the findings of ljenkins et al. I <200lh '). 



Fig. [16] deviations from universality are even greater than simply 
using & = 0.2. The conversion formula is not accurate, nevertheless 
even using the SO halo finder with density threshold A^ir gives 
similar results. From Fig. [16] we can also notice that this type of 
halo detection tends to overcorrect the different mass functions dif- 
ferently for the different cosmologies. The largest deviation occurs 
for L-ACDM, then L-RPCDM and ACDM-W5. This means that 
other effects may contribute to the final shaping of the mass func- 
tion, effects which go beyond the simple spherical collapse model. 
Overall using bvir to account for the virialization process is in- 
correct, for this very reason it will be useful to find an empirical 
relation which can account for it. This will be discussed in the next 
paragraph. 



6.3 An empirical relation: buniv vs. A^ir 

Rather than trying to identify non-linear mechanisms which are not 
well modelled by the spherical collapse we may take a different ap- 
proach and investigate the following question: which virial density 
parameter Auniv or alternatively linking length buniv is needed to 
recover a universal behaviour to numerical precision of the simula- 
tion? 

One pos sibility would be to use a halo mass conversion for- 
mula such as IHu & Kravtsov however this would lead to 
results which are dependent on the specific form of the halo pro- 
file. Therefore we prefer to adopt a more robust approach which, 
although more time consuming, is independent on the halo pro- 
file. We have run a series of 7 FoF halo finders with various linking 
lengths ranging from 0. 15 to 0.21 (in steps of 0.01) for each cosmo- 
logical simulation and redshift. For each measured mass function 
we have computed the deviation from universality and through in- 
terpolation we have determined the linking length parameter which 
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Figure 17. Sensitivity of halo mass functions to the FoF linking length. The 
plot shows the deviations of the mass functions in the plane / — ln{St / a) 
(where 5* = 5c /1. 686) for linking length ranging from 0.15 (lower curve) 
to 0.21 (upper curve) by step of 0.01. This plot is only true for ACDM 
W5 cosmology at z=0. The deviations from universal behaviour also affect 
small masses. 



minimizes the residual with respect to SCDM*. For instance in 
Fig. [it] we plot the deviations measured for the various linking 
lengths in ACDAI-W5 simulations at z=0. Notice that deviations 
are also important at low masses unlike when varying Sc. In this 
specific case the residual is minimal for buniv = 0.187 (between 
the third and fourth curve from the top). For each data-point shown 
in Fig. [15] this procedure allows us to determine the value of buniv, 
or better still {buniv /0.2)~'^ for which a given data-point of each 
mass function lies on the zero residual axis. This provides us with 
an ensemble of {buniv /O-'i)'^ value for each cosmological model 
and at each redshift (2 = 0, 0.25 and 1), for which we calculate 
the average and the standard deviation. So in total we have nine 
estimations that we plot in Fig. [18] as a function of Avir{z)/178 
for each cosmology and redshift. Assuming the conversion for- 
mula A/178 ~ (0.2/6)"' is valid, we can interpret this plot as 
{buniv /O-'i)'^ as a function of (6^i,./0.2)^^ or alternatively as a 
plot of Auniv / AscDM* as a function of Ai,ir/178. Not surpris- 
ing we find that these points are not randomly distributed but are 
compatible with a linear regression. We find 
-3 



0.2 



0.24 X 



178 



+ 0.92. 



(23) 



which we plot in Fig.[T8]as solid line. We can see that a fixed univer- 
sal value of the linking length parameter b = 0.2 (short dashed line) 
does not exist, as it is ruled out at more than 3cr. Similarly the spher- 
ical collapse model prediction Avir/178 ~ {0.2/buniv)'* (dotted 
line) is ruled out, this was expected since the spherical collapse 
is a too simplistic model. Nevertheless even in the context of the 
spherical collapse, accounting for the virialization provide a good 
qualitative description of the dark matter halos. As an example, in 
L-ACDM halos are more "compact" than in SCDM given the dif- 
ferent value of A^i,., hence requiring a smaller linking-length pa- 
rameter buniv as shown in Fig.[T8] The linear regression best-fit lies 
between these two extreme cases. However, this relation should not 
be interpreted as a new form of universality, as indicated by the fact 
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Figure 18. FoF linking length (fe„„i„/0.2)~"^ (or corresponding 
^univ / '^SCDM*) required to explain the deviations from universality 
as a function of virial overdensities from the spherical collapse A„ir/178 
(or corresponding FoF linking length {byir/0.2)~'^) . The symbols are the 
same as for previous graphs. The horizontal dashed line corresponds to the 
usually assumed universality which is strongly ruled out here. The dotted 
upper line corresponds to deviations exactly explained by the spherical col- 
lapse. This is also ruled out. The continuous intermediate line is the linear 
best fit. This correlation indicates that dark energy effects are indeed encap- 
sulated in ^yi^p (and 5c) but in a non trivial way. 



that the slope is neither or 1 . Moreover the dispersion of the points 
around the linear regression is a signature of all the effects which 
contribute to departure from universality such as non- sphericity, 
concentration parameters, halo merging rates, etc... which are not 
modelled by the spherical collapse and which may also depend on 
the cosmological model (i.e. the properties and abundance of dark 
energy). 

Finally, in order to show that Eq. J23b explains most of the 
deviations from a rigorous universal behaviour of the mass func- 
tion to numerical precision, we plot in Fig [T9] the residual of the 
mass function for SCDM*, ACDM-W5, L-ACDM and L-RPCDM 
at z = 0, 0.25 and 1 as a function of 5*/cr , with each mass function 
determined using a linking length parameter buniv obtained from 
Eq. ( 123 1 . The reference cosmology is SCDM* with FoF(fe — 0.2). 
As we can see all deviations reduce to below the 5 percent level, 
which is of order of our numerical precision. 

These results show the importance of taking into account non- 
linear effects in the prescription of the mass function. Again we 
want to stress that we are not claiming to have found a new uni- 
versal behaviour. As already mentioned other non-linear effects, 
besides those modelled by the spherical collapse contributes to the 
mass function. Future works with much smaller numerical errors 
may found departures from Eq. i23l which are correlated with other 
non-linear collapse quantities. Here we have shown that taking into 
account the time of collapse (as encoded by the density threshold 
Sc) and the virialization process (as described by the virial over- 
density Avir) play an important role in determining the halo mass 
function and as such these effects should be included in a theoret- 
ical formulation, also extended to dark energy cosmologies, since 
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Figure 19. Deviations from our reference mass functions (SCDM*, FoF, 
b=0.2) for all cosmologies and redshifts when taking into account the cos- 
mology and redshift dependent time of collapse (by plotting as a function 
of (5» = <5c/ 1.686) and virialization processes (by using a linking length 
b^niv i^vir) given by Eq. l23l where A„ir is given Fig. [5). The convention 
for the cosmologies are the same as in precedent plots: SCDM* (squares), 
ACDM (triangles), L-ACDM (crosses) and L-RPCDM (diamonds). The 
redshifts are z=0, 0.25 and 1 from left to right. The average position of 
the mass functions for all cosmologies and redshifts are well within the 5 
percent deviations band (dotted lines) and are of order of our numerical 
uncertainties. Overall this shows than one can develop a prescription to go 
beyond the universality approximation which accounts for dark energy and 
redshift dependant non-linear effects. 



they greatly improve the agreement between theory and simula- 
tions. 

From a phenomenological point of view, this implies that 
when constraining dark energy models through measurements of 
cluster number counts, rather than using standard universal fitting 
formula of the mass function, it will be preferable to use a prescrip- 
tion which explicitly depends on 5c{z) predicted by a given dark 
energy model. Here, we have shown that Eq. l l22t . with SCDM* 
parameters (given at the beginning of Section |6j provide an accu- 
rate description of the mass function in dark energy cosmologies. In 
addition, one can account for the effect of the virialization by using 
Eq. l l23t . with A„ir value predicted by a given model. In such case 
one has to convert the corresponding mass Mh^^^-^ to the observed 
one usually defined at a given overde nsity Ma ■ These c onver sions 
can be performed using for instance IHu & Kravtsov I l l2003h and 
iLukic etal.1 ( l2009l) . 

A final point concerns whether the non-universality found here 
extends to models with different values. In order to isolate the 
effects of dark energy we have focused on toy-models for which 
we have forced erg to be that of the reference ACDM-W5 cosmol- 
ogy. Would we observe a non-universality also in "realistic mod- 
els" of dark en ergy calibrated on SNIa and CMB data such as those 
considered in l Alimi et al. Il2010l) ? In this case erg differs from one 
cosmology to another, but the deviations from universal behavios 
must be present also for these models. In fact from the analysis 
of the mass function in the WMAP cosmologies presented in Sec- 
tion lS.ll we have shown that erg has very little effect on / (In ct~^). 
Despite having very different ag values the mass function of the 
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WMAP models are universal to numerical precision. This is be- 
cause what really matters is whether cosmological models share 
the same structure formation history or not. Therefore our conclu- 
sions on the non-universality of the mass function can be easily ex- 
tended to cosmologies with different erg. For instance, from Fig. [8] 
and Fig.[TT] it is straighforward to see that the deviations between 
L-RPCDM and ACDM-W5 (we share the same erg) still hold if 
we had confronted L-RPCDM to ACDM-Wl (for which erg is very 
different). Furthermore, having shown that the non-universality of 
the mass function results from the cosmological and redshift de- 
pendence of the past structure formation history, it implies that our 
conclusions do not simply restrict to the models considered here, 
but have more general validity. 



7 CONCLUSION 

In this work we have studied the cosmology and redshift depen- 
dence of the halo mass function through high resolution cosmo- 
logical N-body simulations. By comparing the results of ACDM- 
WMAP calibrated models against toy-models characterized by the 
same distribution of linear density fluctuations today, but with dif- 
ferent expansion histories and growth of the linear perturbations, 
we have been able to infer a number of new and important results. 
Previous works have focused on the universality of the halo mass 
functions and its deviations as a function of redshifts. We have used 
a FoF(6 = 0.2) halo finder to construct catalogues of halos from 
our simulations, and we have limited our analysis to mass ranges 
for which the Poisson noise is below 10% level and halos contain 
at least 350 particles. By focusing on mass function residuals we 
have further limited the effect of numerical systematics uncertain- 
ties to better than 5%. Using such an approach we have been able 
to clearly show that universality does not exists in absolute terms, 
rather it can be verified to numerical precision at 2 = 0, but only 
for those cosmologies which share very similar evolution histories 
at the background and linear level. This is indeed th e case of the 
ACDM-WMAP models, thus confirming past results l lTinker et al. I 
l2008l : ICrocce et alj|2010l) . In contrast our toy-models show depar- 
tures from universality above numerical uncertainties, larger than 
10% at z = 0. 

Using the spherical collapse model as guiding tool, we have 
been able to incontrovertibly identify the non-linear mechanisms 
responsible for such deviations. Firstly, the spherical collapse 
threshold which differs from one cosmological model to another 
is responsible for deviations in the high mass end of the mass 
function. Models with values of 5c lower than the SCDM* pre- 
diction form structures earlier, thus leading to different imprint on 
the present halo mass function consistently with the exponential 
cut-off expected in the Press-Schechter formalism. Indeed account- 
ing for the collapse threshold reduces the discrepancy between the 
mass function of the different models. We have provided a fitting 
formula, based on the Sheth-Tormen functional form calibrated on 
the ACDM-W5 model for which model dependent deviations from 
universality are of ~ 10%, with an explicit dependence on 5c, thus 
it can be efficiently used for a robust model parameter inference 
against halo mass function measurements. 

Nevertheless the residuals lie still above numerical errors and 
are not correlated with the value of 5c predicted by each model, 
suggesting that at least one additional non-linear mechanism is re- 
sponsible for departures from universality. We find that the residual 
deviations at z < 1 are correlated with the redshift evolution of the 
virial overdensity A„ir specific to each model. On the other hand. 



the mass functions tend to that of the standard cold dark matter 
scenario at higher redshifts, which is consistent with the fact that 
all cosmologies tends to a matter dominated universe at z > 1. 
We find an empirical relation between the linking length parameter 
of the FoF algorithm necessary to recover a universal form of the 
mass function and the virial overdensity at a given redshift for each 
cosmology. Such a relation lies between the prediction of purely 
spherical collapse and a fully universal behaviour. It also suggests 
that other non-linear mechanisms probably exists within the nu- 
merical precision of our simulations which contributes to further 
shaping the mass function. 

Furthermore the empirical relation we have found between 
buniv or equivalently A„„i„ and Aui^ may have important im- 
plications for the definition of the physical "frontier" of halos. As 
we have seen here, most of the works in the literature which con- 
cern the mass function use a mass definition based on a constant 
detection parameter (either A or 6). In contrast internal halo pro- 
files are usually defined as a function of A.vir- Henceforth, using 
a value of A„„i„ (or buniv) corresponding to that of the assumed 
cosmology may provide a consistent halo definition applicable to 
both mass function and halo profile. Overall these results point 
to the fact that the collapse threshold and the virialization process 
play an important role in determining the halo mass function, en- 
coding cosmological (dark energy) dependent feat ures which are 
neglected in the standard universal fitting formula jlenkins et al. I 
l200lHWarren et al. II2006I) . and which we have shown to be of lim- 
ited precisions. 
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Alimi et ^ 120101 ) and the halo profile jWechsler et al. 

Dolagetal. 2004"), and which are a direct consequence of the fact 
that the non-linear structure formation conserves a record of the lin- 
ear phase of collapse. Henceforth the halo mass function contains 
specific cosmology dependent features which can be tested through 
observations, provided that predictions from numerical simulations 
or semi-analytical pr escriptions are s ufficiently accurate to account 
for such imprints (see iWu et al.l2O10l) . The results of this work pro- 
vide also an understanding of the deviations from universality that 
are present in realistic dark energy models ( calibrated on CMB and 
SNIa data) as those discussed in lAlimi et al. (2010), and which 
will be presented in an upcoming paper. Finally, we can speculate 
that our findings are relevant also in the context of non-minimally 
coupled inhomogeneous dark energy models, for which the linear 
growth history is scale dependent, thus deviating from that of stan- 
dard LCDM cosmologies. In such a case imprints on the halo mass 
function should be larger, thus needing accurate studies of the non- 
linear structure formation of these scenarios beyond those already 
discussed in the literature. 
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